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Abstract 

The spectrum of the T system is investigated using the Nonrelativistic Lattice QCD 
approach to heavy quarks and ignoring light quark vacuum polarization. We find good 
agreement with experiment for the T, T', T" and for the center of mass and fine struc- 
ture of the Xb-states. The lattice calculations predict bb D-states with center of mass 
at (10.20 ± 0.07 lb 0.03)GeV. Fitting procedures aimed at extracting both ground and 
excited state energies are developed. We calculate a nonperturbative dispersion mass for 
the T(IS') and compare with tadpole-improved lattice perturbation theory. 

PACS numbers: 12.38. Gc, 14.40. Gx, 14.65.Fy, 12.39.Hg 



1 Introduction 

Hadrons containing one or more heavy quarks have been the focus of intense investigations by 
lattice gauge theorists in recent years. One motivating factor is that these systems are also 
being studied extensively by experimentalists trying to nail down the remaining parameters in 
the Standard Model. Nonperturbative QCD results are needed in many instances, to convert 
experimental numbers into determinations of fundamental parameters or to test the Standard 
Model. The lattice approach to nonperturbative QCD is now starting to yield reliable numbers 
for several of these crucial inputs. Part of the activity has been in heavy-light systems, focusing 
on leptonic and semi-leptonic decays of heavy-light mesons (the B's and D's) and on neutral 
meson mixing |T| . Another area of investigation, which is also the focus of the present article. 
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has concentrated on heavy-heavy systems such as the J/\l/ and T famihes. Studies of the latter 
systems have aheady lead to the most accurate lattice determinations of the strong coupling 
constant, 0, and of the b-quark pole mass Mf, [^. In heavy-heavy systems one can 

take advantage of the fact that only heavy quark propagators are required to do high statistics 
simulations at only modest computational cost (of course only once the gauge configurations 
have been created). This coupled with the wealth of experimental data on quarkonium allows 
one to carry out stringent tests of computational methods employed by lattice gauge theorists. 
These systems may also be the place where effects of quenching can be studied quantitatively. 
We report here on a study of the T system using the Nonrelativistic QCD (NRQCD) p, ^ 
approach to heavy fermions. Our goal is to start from a first principles QCD hamiltonian and 
show that it can reproduce the T spectrum. Along the way we develop and refine methods 
to analyze numerical data, methods that we hope will be useful in other lattice calculations 
as well. For instance, we find that extracting excited state energies is straightforward using 
our fitting procedures. Our simulations also serve to test perturbation theory on the lattice 
in a new setting, through comparisons of nonperturbative simulation results for the T kinetic 
mass with perturbative formulas. The investigations in this article provide the foundations for 
our determination of the b-quark pole mass. The calculations are described in a separate 
publication]^. We emphasize that NRQCD provides an extremely efficient way to obtain 
realistic and accurate heavy quark propagators. The prospect of NRQCD having impact not 
only in investigations of heavy-heavy, but also of heavy- light systems looks very promising [^]. 

The b-quarks in the T system are quite nonrelativistic, with v"^ ~ 0.1. The splittings 
between spin averaged levels are around ~500 MeV ((^(Mftf^)), which is much smaller than 
the mass ((9(2Mf,)), indicating that a systematic expansion of the QCD hamiltonian in powers 
of v"^ is very appropriate here. The continuum action density, correct through 0{MbV^), is 
given by 

with 



TTCOnt. 

" ~ 2A.f» 



(2) 



if), i/j^ are two component Pauli spinors and at tree-level we have Cj = 1 for all i. Previous 
NRQCD studies ^ have used the leading order Hamiltonian {6H = 0) or the leading order 
plus the <T ■ B term and gave encouraging results. Here we include all the 0(Mbf'^) terms. 
This means systematic errors due to relativity will be of 0{MbV^) ~ 5MeV, which is about 
1% of a typical radial or orbital excitation energy and 10% of a typical spin sphtting. 
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In our lattice simulations there will be other sources of systematic errors. Finite volume 
errors are not a problem here since the T's are smaller than regular light hadrons. Finite 
lattice spacing errors in the fermionic action can be corrected for order by order, similar to 
the systematic f ^ corrections. So the main remaining sources of systematic errors come from 
the gauge configurations, namely 0(a.^) errors from using the standard Wilson gauge action 
and quenching errors due to the absence of light quark vacuum polarization. From potential 
model calculations we estimate the latter errors to be the dominant ones. One consequence of 
quenching will be that inverse lattice spacings, a^^, extracted from different observables will 
not agree with each other in general and one will have to make some choice when presenting 
dimensionful results. We find that our two basic splittings, the IS- IP and the 1S-2S T energy 
level splittings, give a~^'s that differ by one to two a. We use an "average" lattice spacing in 
our dimensionful plots. 

Finally we need to discuss the number of parameters in the NRQCD action. In addition 
to the bare mass, M°, and the gauge coupling, g, one has the q's. We work with the q's 
set to their tree-level values q = 1 while at the same time "tadpole-improving" the lattice 



version of the NRQCD action |T^. This ensures an optimal perturbative scheme so that one 
can expect renormalization effects to be small. We find that in practice, tree-level values 
give the correct P-state fine structure splittings, giving us confidence that setting q = 1 
is a viable approach. The final parameter we must consider is the zero of energy which 
is used to relate NRQCD energies with absolute, relativistic energies (this term is usually 
omitted from NRQCD actions). As with the Cj's, it can be fixed at tree- level to equal M°, or 
calculated in perturbation theory. It can also be calculated non-perturbatively by requiring 
that the dispersion relation of the T be Lorentz invariant, up to the order in at which 
we are working. We find excellent agreement between perturbative and non-perturbative 
determinations, further encouraging us that perturbation theory is working. We stress that 
the only free parameters which we tune to match experimental results are those appearing in 
the original QCD action, M° and g. In other words, this is a first-principles QCD calculation, 
not a QCD inspired phenomenological model. The coupling g is eliminated as a free parameter, 
in the usual way, when we fix the scale = a~^{g) to match the IS - IP and/or the IS 
- 2S splitting. M° is tuned so that the simulated kinetic mass for the T agrees with the 
experimental T(IS') mass. 

Our results for the spectrum of the T system are shown in Fig. |1] and Fig. 0. We use 
= 2AGeV, which is an average between the inverse lattice spacings obtained from the T 
1S-2S and IS-IP splittings (the error in this estimate for is at the 4% level, details are 
given in section 4). One sees that the general features of the known spectrum are reproduced 
nicely 0]. Fig. |l] shows the T, rjb and singlet P- and D- states and Fig. |^ shows the P-state 



fine structure. In both figures the errors reflect statistical errors plus some systematic fitting 
errors. We do not show estimates of systematic errors due to quenching or the effects of 
uncertainties in the scale a~^. 

Both the rib and the D-state (center of mass) are predictions of the theory. We find an 
r]b state at (9.431 ± 0.005 ± 0.001)GeV and D-states with center of mass at (10.20 ± 0.07 
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Figure 1: NRQCD simulation results for the spectrum of the T system including radial exci- 
tations. Experimental values (dashed lines) are indicated for the triplet S'-states, and for the 
spin-average of the triplet P-states. The energy zero from simulation results is adjusted to 
give the correct mass to the T( 1^5*1 ). 

± 0.03)GeV. These numbers include the dominant statistical and/or systematic errors other 
than those due to quenching. For the D-states the first error corresponds to the statistical 
error in fitting the D mass in lattice units. For the rjb state this error is neghgible and so 
the first error quoted there is the systematic error from neglected higher order relativistic 
terms and finite lattice spacing corrections. In both cases the second error arises from the 
uncertainty in the value of a~^. We expect the S-states and hence also the T - rjh splitting to 
have noticeable quenching errors. Spectrum calculations with dynamical gauge configurations 
are already underway. It will be interesting to compare the quenched and unquenched spectra. 

In the rest of the article we give more details of our analyses, starting with the quark 
propagator calculations and the meson correlations in the next section. Section 3 includes a 
lengthy explanation of our fitting procedures. The main message there is that, instead of going 
to large timeslices in search of a plateau, we have worked with high statistics correlations on 
shorter lattices (24 time slices) and carried out simultaneous multi-exponential fits to several 
correlations at a time. We also devised methods to extract splittings directly. Our statistical 
errors for spin averaged level splittings are about 20 - 30 MeV, for P-state fine structure 
splittings about 5-10 McV and for the S-state spin splitting down to O.SMeV. In section 4 
we discuss determination of a^^ and comparison with experiment. We also present results for 
meson wave functions at the origin. Section 5 describes the extraction of a "kinetic mass" , 
aMkin, for the T(15') from correlations with momenta and comparisons with perturbation 
theory. Section 6 gives a brief summary. 
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Figure 2: Simulation results for the spin structure of the lowest lying P-wave states in the 
T family. The dashed lines are the experimental values for the triplet states. Energies are 
measured relative to the center of mass of the triplet states. 

2 The Simulation 



2.1 Quark Propagators 

Quark propagators in lattice NRQCD are determined, in a single pass through the gauge-field 
configuration, from evolution equations that specify the propagator for t > in terms of its 
value at t = 0. Various evolution equations have been suggested in the past. Currently we 
use the equation proposed in 0, modified slightly for improved efficiency. Our propagators 
are defined by the equation 



1- 



2n 



Ua 1 



aHr ^ 



2n 



(3) 



where = for t < 0. For numerical work it is convenient to rewrite this equation in the 
form: 



Gi 
Gt+i 



aHi 



2n 
aHo 



Ulil 



2n 



x,0 



2n J ^ V 2n 
On the lattice, the kinetic energy operator is 



l-a5H)Gt (t>l). 



2Mfe°' 



(4) 



(5) 
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and the correction terms are 



-Cs^AttttO- ■ (A X E - E X A) - Ci-^ a ■ B 

_ . - ,n - Cfi . „ (6) 



The last two terms in 5H come from finite lattice spacing corrections to the lattice laplacian 
and the lattice time derivative respectively. A is the symmetric lattice derivative and A*^^) is a 
lattice version of the continuum operator Df. We used the standard cloverleaf operators for 
the chromo-electric and magnetic fields, E and B. The parameter n is introduced to remove 
instabilities in the heavy quark propagator caused by the highest momentum modes of the 
theory. For our simulations at /3 = 6.0 and with bare masses relevant for the T system, we 
set n = 2. 

As mentioned in the Introduction, we tadpole-improve our lattice action by dividing all 
the U's that appear in E, B, and the covariant lattice derivatives fields by uq, the fourth root 
of the plaquette||lO|. This is most easily done as the t/^'s are read by the simulation code. 



The effect is to transform the operators that appear in the evolution equations as follows: 

ul Ul/uo (7) 

E E/iuo)' (8) 

B B/M^ (9) 

A^,G{x) — y (U^{x)G{x + fi)-Ul{x-fi)G{x-fi))/2uo (10) 

A.fG{x) (u^{x)G{x + fi)+Ul{x-fi)G{x-fi))/uo-2G{x) (11) 

with 



A(^) = Y.^f (12) 
AW = i:{Aff (13) 

Tadpole-improvement of the action allows us to work with tree-level values for the Cj's in 5H 
without having to worry about large renormalizations. Hence our lattice action depends only 
on two parameters, the bare mass M° and the QCD coupling constant g. We have collected 
data for three values of the bare mass, aM^ = 1.71, 1.8 and 2.0, all at /5 = 6.0. For this P, 
uo = 0.878. 

We computed our quark propagators for aM° = 1.8, 2.0 using a different evolution equa- 
tion: 

U,Gt+r-[l--^] {l-a5H)[l-—^] Gt = 6,,oSt,o- (14) 



2n J ' ' \ 2n 

Propagators from this evolution equation are the same as those from the other (Eq. (|^)), but 
with extra factors of (1 — aHo/2n)~^ in both the source and sink. These extra factors cannot 
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affect the meson spectrum, but they do modify the quark's wavefunction normahzation in 
order f^. For this reason our other formulation (Eq. (^)) is superior; with it, NRQCD quark 
fields are normalized to unity up to corrections of order and w^. 



2.2 Meson Correlation Functions 

Once one has the quark propagators it is straightforward to obtain meson propagators. Let ip'^ 
and denote fields that create a heavy quark or heavy anti-quark respectively. The following 
interpolating operator creates a meson of momentum p. 

^H^i) r(fi - X2) xHx2) e'i-(^i+^2) (15) 

Xl,X2 

where the "meson operator" = r(a;i — X2) = fl (pdxi — X2\) ■ The operator is a 2 x 2 matrix 
in spin space and generally includes derivatives acting on the radial function (j){r). Using 
translation invariance, we eliminate the summation over the initial antiquark position. The 
meson propagator is then 



=Y.Tr [Glim) r(^'=)^(fi - y2) Gt{yi] 
with 

Gt{y)^Y.Gt{y-S)V^''\S)^'^' (17) 



and the trace is over spin and color. In the above equations we distinguish between F'-*'^-' and 
F("^), i.e. the smearing at the source or sink. Gt{x) can be obtained directly using Eq. (1) 
with ^ T^^^\x) e*2'^. In the future we will often refer to the smeared propagator G as the 
quark propagator. The convolution in Eq. (|T6|) is evaluated using fast fourier transforms. 



In Table |l| we list the zero momentum bb meson states studied in the current project 
together with their corresponding "meson operators", F(a;). We choose 4>nsc{^) for ^sc = 1, 2, 3 
to correspond to Richardson potential radial wave functions for the S-, P- or D-state ground 
and excited levels. We also used 5-function local sources and sinks. These are referred to as 
^sc = loc in Table |I[ For each set of meson quantum numbers we evaluated correlations with 
all possible independent smearings,(nsc, ns^) = (/oc, /oc), (/oc, 1), (/oc, 2)... etc., at the source 
and the sink. So for the S-states we obtained a 4 x 4 matrix of correlation functions and 
for the P- and D-states 3x3 matrices. In addition we looked at S-state mesons with small 
momenta. For those we used (pioc and 0i. 

Table |l| also lists the continuum quantum numbers J^'^ and the lattice cubic group repre- 
sentations for our meson states. This allows us to see which states are expected to mix with 
each other. Mixing will occur both due to relativistic corrections and due to lattice artifacts. 
For instance, since L is no longer a good quantum number in a relativistic theory the ^Si 
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Table 1: Meson Operators. Aj denotes the symmetric lattice derivative and Dij = AjAj — 
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and the ^Di states will mix, both being J^'~"' = 1 . This mixing, which happens even in the 
continuum limit, is suppressed by f^. On the lattice one can also have mixing between dif- 
ferent J-states that fall into the same lattice representations. Examples are mixing between 
^D2,T2 and ^-02^2 or between ^D^Ti, ^DiTi and ^5*1. We have measured cross-correlations 
between these states, but postpone their analysis for future work, concentrating here on the 
spin-averaged D-states. 

In the expression Eq. (|15[) one has sums over color and spin degrees of freedom. One 
could calculate quark (and antiquark) propagators separately for each color and spin quantum 
number at the source. We have done so for the color degrees of freedom and verify a reduction 
of statistical errors by v^S compared to when only one value for the initial quark and antiquark 
color was used. As far as spin is concerned we decided to save on CPU time by setting the 
initial quark and antiquark spins equal to 1. This means that at the source we are sensitive 
only to the 1-1 component of the 2x2 spin matrix in V'^^^'^ and mesons of definite quantum 
numbers are projected out at the sink. From Table one sees that groups of mesons such as 
^5*0 and "^Siz , or ^Pix , ^Piy and ^P2T2Zx etc. have the same 1-1 component of T^^^^ up to 
normalization. For each group, the meson correlations for its members can be obtained from 
one common quark propagator (this must be repeated for each smearing function, at the 
source), and are highly correlated. We have taken advantage of these strong correlations to 
reduce statistical errors in our fits for hyperfine and fine structure splittings. We worked with 
13 zero momentum quark propagators with S- or P-state smearing at the source and four S- 
state propagators with momentum. Out of these quark propagators 129 S- and P-state meson 
correlations were evaluated using different r'^'^*^-' and many combinations of smearing functions 
at source and sink. For aM° = 1.71 we also evaluated an additional 15 quark propagators 
with D-state smearing. To date we have only analyzed combinations of these D-state quark 
propagators giving rise to the mesons. With 5 polarizations and 9 different source-sink 
smearing combinations, this means a total of 45 D-state meson correlations. 



3 Data Analysis and Fitting Results 

Our calculations used (3 = 6.0 quenched gauge field configurations on 16^ x 24 lattices provided 
by Greg Kilcup and his collaborators. We worked with an ensemble of 105 independent 
configurations that were gauge-fixed to Coulomb gauge. For each configuration we selected 8 
different origins on the first time slice for our quark (antiquark) propagators. Simple tests, 
in which data from different origins but from the same configuration were binned, indicated 
that different origins lead to independent propagators. For our aM° = 1.71 P-states we also 
ran with 8 origins on timeslice 12. This data was binned together with the timeslice 1 data 
for each configuration and spatial origin. Hence for each of the 129 (or 129 + 45) meson 
correlations, discussed above, we worked with 105 x 8 = 840 measurements. In Fig. 3 and 
Fig.4 we show some examples of effective mass plots for our data. The errors are bootstrap 
errors. Fig. 3 shows aM° = 1.8 data for ^5*0 states. The effective mass plots are arranged as a 
3x3 matrix corresponding to the 9 source-sink combinations {ugc, ngk) with n^c, = loc, 1,2. 
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Figure 3: ^5*0 Effective masses by (source, sink). 



Fig. 4 shows similar plots for aM^ = 1.71 ^Pi data. We have averaged over ^Pix, ^Piy and 
^Piz. This hb state has not been observed yet experimentally but there are strong theoretical 
reasons for believing that it lies close to the center of mass of the ^P levels. Hence we will 
sometimes refer to the ^Pi level as the "spin averaged" P-state. One sees from Fig. 3 and Fig.4 
that our S-statc correlations with Ugk = loc and = 1 sinks have excellent statistics. The 
S-state correlations with excited state smearing Ugk = 2,3 and the P-state data have good 
to reasonable statistics. In the effective mass plots, truncations at large t means that signal 
to noise in the original data was worse than 3 : 1 beyond that point or that the correlation 
had switched sign. We have used the naive definition m^ff{t) = log{G{t)/G{t + 1)), although 
a more sophisticated version could be used when one deals with off-diagonal correlations in 
which some amplitudes can come in with negative signs. Our plots provide a rough assessment 
of the quality of the data; we do not use them in our fits. 
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3.1 Fitting Procedures for ^5*1 and Singlet P- and D-States 

We investigated a variety of fitting procedures to extract the spectrum from the meson corre- 
lations. For the S-states and singlet P- and D-states we used simultaneous multi-exponential 
fits to several correlations in order to obtain the ground state and one or two excited state 
energies. We employed mainly two such fitting procedures and got consistent results from 
both of them within one a. In the first procedure, we fit simultaneously to a matrix of correla- 
tions with Use-, nsk = 1; 2, 3 for S-states and n^c, Ugk = 1,2 for the spin singlet P- and D-states. 



Similar fitting methods have been discussed by other groups in the past[0. In the matrix fits 
each meson correlation is fit to the expression 

^ exp 

Gmeson{nsc,nsk]t) = a^Usc, k) a* {usk , k) e~^'''\ (18) 
fc=i 

Using charge conjugation symmetry one can argue that, for equal mass quarks, the coefficients 
a{j, k) can be chosen to be real numbers. From the effective mass plots of Fig. 3 and Fig.4 it is 
clear that if one wants to make use of most of the data one does not want to move too far out 
with tmin in search of a plateau. One wants instead to fit to several exponentials working with 
as small a tmin value as possible while still maintaining good chi-squareds. Multi-exponential 
fits to a single correlation, however are usually tricky and unstable. It is much easier to 
do multi-exponential fits to several correlations simultanously, especially if one makes sure 
that each exponential has a significant amplitude in at least one correlation. Although the 
number of parameters in our fits was sometimes large (we were typically doing 6-, 9-, 12- and 
16-parameter fits) they were highly constrained and we did not run into stability problems. 
Since we had a large number, 840, of measurements to play with we also did not have to worry 
about the size of the covariance matrix becoming too large. We used Ngr^p = 3 for the 3x3 
matrix fits and Ng^p = 2 for the 2x2 matrix fits. 

The second procedure for doing mult i- correlation multi-exponential fits was to take a set 
of smeared-local correlations, {nsc,nsk) = (nscjoc) with Ugc = 1,2,3 or n^c = 1,2 and fit 
them simultaneously to the same set of energies, Ek- These smeared-local correlations have 
the smallest statistical errors. They are fit to 

Nexp 

GmesoninscJoC;t) = J2 K^^sc, k) e~'^'' \ (19) 
k=l 

Note that this is the same ansatz as Eq. (|l^), with b{nsc,k) = a{nsc, k)a*{loc, k). We used 
Nexp = 2 or 3 when fitting two correlations and Nexp = 3 or 4 for three correlations. We 
use different sets of correlations for fitting procedures 1 and 2. So the fact that, as we will 
see, consistent energies are obtained from the two procedures gives us confidence in the final 
results. 

In fitting correlations the delicate question is always how to pick the range tmin/tmax in i 
over which to fit. Needless to say, in our multi-exponential fits the larger Nexp the smaller we 
can make tmin- On the other hand for fixed N^xp when one increases tmin there could come a 
point when the signal for the higher exponential is no longer above the noise leading to larger 
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Table 2: Examples of simultaneous multi-exponential fits to two and three ^5*1 smeared-local 
correlations. 

errors in the energies of even the lower lying levels. To illustrate some of these issues we show 
in Table examples of fits to aM° = 1.71 '^S'l-states using procedure 2. We tabulate results for 
different number of exponentials and correlations and show the dependence on t min /t max- The 
last column gives Q, the "goodness of fit", i.e. the probability that fluctuations in correctly 
modelled data will generate a greater than that of the fit. One usually desires Q > 0.1, but 
slightly smaller values do not necessarily rule out a fit. Table ^ starts with two exponential 
fits to two correlations with (nsc,nsk) = (l,loc) and (2,loc). One sees that for tmin of 6 to 8 
good fits are obtained to the ground- and first excited level. For tmin = 4 the Q- value starts to 
deteriorate. Adding a third exponential allows one to go down to tmin of 2. The third energy, 
E3, however is not reliable yet until we include another correlation with 03 smearing. Such 
fits, i.e. simultaneous fits to three correlations, Gss^i^, ^oc; t), 6^35^(2, loc; t) and G35^(3, /oc; t) 
are also shown in Table 

In Table |^ we give examples of 2 x 2 and 3x3 matrix fits (procedure 1) to ^Si levels 
for several tmin/tmax- The results are consistent with Table ^ but with generally worse Q- 
values. For the 3x3 matrix fits one sees that the signal for the second excited state has 
disappeared once tmin > 6. The errors in Table || and ^ are obtained by the criterion that 
6x^ = 1- We have checked that bootstrap errors agree with these errors to within ±10%. This 
indicates that the statistical fiuctuations in our correlations are close to being Gaussian. In 
Fig. 5 and Fig. 6 we show "effective amplitude" plots (for Gmesonit) ■ e"^^'*) corresponding to 
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-^exp '^min /^max ClEi 


aE2 




Q 


fits to 


2 


5/24 


0.4537(7) 


0.694(7) 




0.10 


(1,1), (1,2) 




5/22 


0.4539(7) 


0.693(7) 




0.17 


(2,1), (2,2) 




6/24 


0.4536(7) 


0.708(9) 




0.24 






6/22 


0.4538(7) 


0.708(9) 




0.41 






6/16 


0.455(1) 


0.71(1) 




0.21 






7/24 


0.4536(7) 


0.70(1) 




0.17 






7/22 


0.4539(7) 


0.70(1) 




0.31 






8/24 


0.4534(7) 


0.71(1) 




0.13 






8/22 


0.4537(8) 


0.71(1) 




0.24 


fits to 


3 


4/24 


0.4540(6) 


0.702(6) 


0.88(2) 


0.04 


(1,1), (1,2) 




5/24 


0.4538(6) 


0.697(8) 


0.86(3) 


0.06 


(2,1), (2,2) 




5/22 


0.4540(7) 


0.696(8) 


0.86(3) 


0.06 


(1,3), (3,1) 




6/24 


0.4536(7) 


0.710(9) 


0.90(7) 


0.14 


(2,3), (3,2) 




6/22 


0.4538(7) 


0.709(9) 


0.89(7) 


0.18 


(3,3) 




6/16 


0.455(1) 


0.711(9) 


0.90(9) 


0.05 






7/24 


0.4537(7) 


0.67(3) 


0.74(2) 


0.22 






7/22 


0.4540(7) 


0.67(3) 


0.74(2) 


0.31 






8/24 


0.4534(7) 


0.64(5) 


0.724(9) 0.49 


Table 3: Examples of 2 x 2 and 3x3 matrix fits to ^5"! correlations. 


Fit 




k a{nsc,sk = 1, k) 




--2,k) 


^{nsc,sk = 3, k) 


Nf,xp — 2 fits to 


6/22 


1 


0.877(5) 


-0.065(3) 






(1,1), (2,2) 




2 


0.18(1) 


0.79(3) 






(1,2), (2,1) 


7/22 


1 


0.878(5) 


-0.066(3) 










2 


0.18(2) 


0.78(3) 






Nexp = 3 fits to 


6/22 


1 


0.877(5) 


-0.065(3) 




-0.015(1) 


(1,1), (1,2), (2,1) 




2 


0.18(1) 


0.79(2) 




0.07(4) 


(1,3), (3,1), (2,2) 




3 


0.094(9) 


-0.1(1) 




0.8(2) 


(3,3) 


5/22 


1 


0.878(4) 


-0.066(3) 




-0.014(1) 






2 


0.17(1) 


0.75(2) 




0.03(4) 






3 


0.09(1) 


0.09(8) 




0.70(6) 



Table 4: Examples of fit results for amplitudes a{nsc,sk, k) 
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Fit 


^min/^max 


k 


b{nsc = 1, k) 


b{nsc = 2, k) 


b{nsc = 3, k) 


Nexp — 2 fits to 


6/24 


1 


0.135(2) 


-0.0102(8) 




(l.loc) and (2,loc) 




2 


0.032(4) 


0.114(2) 




Nexp = 3 fits to 


3/24 


1 


0.135(2) 


-0.0105(8) 




(l,loc) and (2,loc) 




2 


0.028(6) 


0.100(8) 








3 


0.030(9) 


0.045(12) 




Nexp — 3 fits to 


7/24 


1 


0.135(3) 


-0.0104(8) 


-0.0023(4) 


(l,loc) , (2,loc) 




2 


0.02(4) 


0.14(8) 


-0.01(3) 


and (3,loc) 




3 


0.01(6) 


-0.04(9) 


0.11(2) 



Table 5: Examples of fit results for amplitudes b{nsc, k) 





Nexp 


tmin/tmax 


aEi 


aE2 


Q 


2x2 matrix 


2 


8/24 


0.630(4) 


0.82(3) 


0.15 


fits 




7/24 


0.630(3) 


0.81(2) 


0.17 






6/24 


0.632(5) 


0.84(1) 


0.03 






8/14 


0.625(5) 


0.83(1) 


0.19 






7/14 


0.626(4) 


0.80(2) 


0.22 


fits to (l,loc) 


2 


8/24 


0.630(6) 


0.85(3) 


0.008 


and (2,loc) 




7/24 


0.628(5) 


0.82(2) 


0.008 






8/14 


0.628(10) 


0.82(3) 


0.23 






7/14 


0.626(8) 


0.81(2) 


0.35 




3 


7/14 


0.629(16) 


0.88(19) 


0.22 






4/14 


0.626(9) 


0.79(6) 


0.40 



Table 6: Examples of fits to ^Pi correlations 
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the Nexp = 3, tmin/tmax = 5/24 fit to two Correlations in Table ^ and the 3x3 matrix fit with 
train/ t-max = 7/22 in Table ^. One sees that most of the correlations are reproduced reasonably 
well by our fits. 

Information on how well our smearing functions are doing is contained in the fitted ampli- 
tudes, the a{n, k)'s and b{n, A;)'s of Eq. (p!8[) and Eq. ([T9|). These measure the overlap between 
our smearing functions and energy eigenstate wave functions. In Table ^ we show a{n, k) 
for several matrix fits of Table ^ and similarly in Table ^ we show b{n, k) of smeared-local 
fits. One sees for instance that with nsc,sk = 1 , the amplitude for the ground state (k = 1) 
dominates, whereas for nsc,sk = 2 or Usc^sk = 3 one has the largest overlap with, respectively 
the first or second excited state (k=2 or k=3). Our Richardson potential smearing functions 
are doing a reasonable job in focusing correlations onto the right levels. One might do better 
by using our own simulations to provide good smearing functions. We are planning to do this 
in future simulations. 

The multi-correlation multi-exponential fits described above were used to obtain energies 
for the ?7b, T, the hf, levels and the singlet D-states. The r]b fits are similar to those shown for 
T states. The hfy and the D-state correlations are more noisy. We give examples of fits for the 
hb levels in Table ^ using 2x2 matrix fits (with ijisc^nsk) = (1,1), (1,2), (2,1) and (2,2)) and 
using simultaneous fits to {risc^nsk) = (l,loc) and (2,loc). We will give summaries of all our 
fits for different M° below, where we also convert dimensionless numbers into real energies in 
GeVs. Here we go on to describe fitting procedures used to obtain fine-structure splittings 
between the Xb states, the T — rjj, splitting and the sphtting between states with zero and 
nonzero momenta. 



3.2 Fitting Procedures for Spin Splittings 

One test of the NRQCD effective action is to see how well the fine and hyperfine structure 
in the bb system can be reproduced. These splittings will be the ones most sensitive to the 
coefficients Cj in the action and the question is whether, with tadpole-improvement, tree-level 
values of Cj = 1 are adequate. The splittings between the Xb states are a few tens of MeV. 
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0.86 
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0.83 
0.82 
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Figure 6: ''iS'i Effective amplitudes G{t) ■ e^^ * for tliree by three matrix fits witli 
Off-diagonal fits are averaged over source and sink. 



^max — 22. 
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From Table § one sees that direct determination of P-state levels have errors that are at best 
0.005 - 0.010 in dimensionless units. With an inverse lattice spacing of ~ 2.4 GeV (see 
next section) this corresponds to an error of 12 to 24MeV and it seems marginal whether 
we would be able to resolve fine structure if each P-level were fit independently. To get 
around this problem one can take advantage of the fact that correlation functions for different 
mesons on the same configuration can be highly correlated. This allows for a direct fit to 
the mass splitting between them. We have employed two methods for extracting splittings, 
6E = E^- , from two very correlated propagators for mesons, A and B, having different 
quantum numbers. We call them the "ratio" method and the "correlated- 6E" method. 

In the ratio method we take our 840 measurements for each meson A and B and create a 
jackknife ensemble of ratios of correlations. This is then fit to a single exponential to determine 
6E. In the correlated- 6E method we pick a set of correlations for each meson and fit the two 
sets simultaneously using 

k=l 

c^K„ 1) e-(^i*+^^)-* + 5 cbK„ k) e-^^* (20) 

k=2 

Of the two methods the second is more general, since it is straightforward to include higher 
exponentials and handle spin splittings of excited states. We have estimated the effects of 
higher states on the ratio method by fitting the jackknifed ratio ensemble to 

Rattoit) = A^e~^^ ' ^ ^ (21) 

We did not find a signal for the correction terms and the value for 6E was consistent with 
those from naive one-exponential fits only with slightly larger errors. We give examples of fits 
to the T - rib splitting and one of the Xb splittings in Tables |^ and |^. For the ratio method we 
use (l,loc) correlations and for the correlated-5i? method (l,loc) and (2,loc) correlations for 
each meson. We see that it is possible to get T — r/b splittings with 0.5MeV and Xb splittings 
with 5 - lOMeV statistical errors (again using of 2.4GeV). 



3.3 Fit Results 

Using the fitting procedures of the previous subsections we obtained energies for the T, T', T" 
levels, the hi,, h[ levels, the T — rib splitting and for the splittings between the xt levels. For 
aMb = 1.71 we also have results for D-states (to date we have only looked at averages over 
the five polarizations of the state and leave more detailed studies of D-state fine structure 
and mixing with S-states to future work). Our estimates for the energies and splittings in 
dimensionless units are shown in Table ^ for several M° values. For comparison we also 
show results for data without relativistic and finite lattice spacing corrections {6H = 0) and 



G meson AiP'sci loC, t) 
G meson b{J^sci loC, t) 
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N 


tmin / tmax 


a5E 


Q 


Ratio Method 


1 


18/24 


0.0123(2) 


0.21 




1 


16/24 


0.0123(2) 


0.21 




1 


12/24 


0.0126(2) 


0.09 




1 


8/24 


0.0129(1) 


0.005 


correlated 


2 


8/24 


0.0123(2) 


0.39 


SE fits 


2 


7/24 


0.0123(2) 


0.49 




2 


6/24 


0.0122(2) 


0.48 



Table 7: ^Sq - ^Si splitting 





N 

^ ^ exp 


tmin/ tmax 


a5E 


Q 


Ratio Method 


1 


10/24 


0.021(4) 


0.77 




1 


9/24 


0.020(4) 


0.79 




1 


8/24 


0.017(3) 


0.70 




1 


7/24 


0.016(2) 


0.76 




1 


6/24 


0.017(2) 


0.77 


correlated 


2 


8/24 


0.023(7) 


0.52 


SE fits 


2 


7/24 


0.021(6) 


0.55 




2 


6/24 


0.018(5) 


0.64 



Table 8: Example of P-state Fine Structure: — ^Pq sphtting 
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aM^ = 1.71 


1 

1.8 


2.0 


1.8 [oH = 0) 


1.8 (mo = 1) 


rSi 


0.4534(8) 


0.4505(10) 


/~\ A A A / 'I \ 

0.444(1) 


0.448(1) 


1.097(1) 




0.695(10) 


0.69(1) 


0.68(1) 


0.69(1) 


1.335(8) 


3-^Si 


0.82(5) 


0.83(4) 


0.83(5) 


0.84(5) 


1.48(3) 


1 1 D 


O.o2o(o ) 


0.62 /[iZ ) 


0.6i9(i0 ) 


0.6o4(io j 


1.27(1) 


2iPi 


0.81(3) 


0.80(4) 


0.79(4) 


0.81(4) 


1.43(3) 




0.76(3) 










Splittings: 












3 c 1 c 
Oi — Do 


0.0123(2) 


0.0116(2) 


0.0106(2) 




0.0049(1) 


'P2 - 'Po 


0.020(4) 


0.021(4) 


0.018(4) 




0.011(3) 


'P2 - 'Pi 


0.008(2) 


0.008(3) 


0.008(3) 




0.0047(14) 



Table 9: Fit results for dimensionless energies and splittings, aE and a6E. 

for another set of data without tadpole- improvement {uq = 1). The errors in Table |^ differ 
slightly from the purely statistical errors of previous tables. In those tables one saw that the 
central values and statistical errors depend sometimes on the fitting procedure used (and on 
tmin/tmax) and wc have tried to take that into account in Table ^ Only fits with good Q- values 
( Q > 0.2 ) and effective amplitude plots were included in these considerations. The reader 
may want to compare earlier tables with the aMj^ = 1.71 column in Table ^. 

There are several features worth noting in Table |^. For instance, the splittings between 
1^5*1 and l^Pi (spin averaged P-state) as well as between 2^5"! and l^S*! levels are insensitive 
to M° within errors. This is known to hold in the real world when one compares splittings 
in the J/"^ and T systems. Also, for spin averaged quantities the 6H = results are almost 
indistinguishable from results with higher order corrections, except maybe for a slight lowering 
of the IS-level. Considering next the spin dependent splittings, one sees that the ^^i-^^o 
splitting does depend on M°, as expected. More striking is the sensitivity of spin splittings 
to tadpole- improvement. The splittings are reduced by a factor of ~ | without tadpole- 
improvement (i.e. if Mo = 1). As we shall see, the tadpole- improved results agree well with 
experiment. Tadpole- improvement of the lattice action appears crucial if one wants to work 
with tree-level values q = 1. Otherwise spin splittings are badly underestimated. 



4 Comparisons with Experiment 

To compare our results with experiment we must convert from dimensionless lattice units to 
physical units by setting the scale a~^. Before doing so, it is worthwhile reminding ourselves 
of systematic errors still contained in our simulations. The largest source of systematic er- 
rors, we believe, comes from the quenched approximation. One expects, among other things, 
corrections to both the IS-IP and the 1S-2S splittings, and inverse lattice spacings obtained 
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by fitting these splittings to experimental data should differ from each other. For instance, a 
bootstrap estimate for the ratio of splittings, {E{2^Si) - E{l^Si))/{E{l^Pi) - E{l^Si)) = R, 
gives R = 1.41(7) for aM^ = 1.71; R = 1.38(8) for aM^ = 1.8; and R = 1.35(8) for aM^ = 2.0. 
The experimental value for this ratio is Rexp = 1-28 and one sees a possible 1 to 2 a deviation 
in our simulations. Inverse lattice spacings obtained from matching the S-P splitting to the 
Particle Data book will differ from that obtained using the 1S-2S splittings by 1 or 2 a. Taking 
splittings from Table |^ that is indeed what we find. Work is underway to repeat simulations of 
the T system with unquenched gauge configurations. It will be interesting to see how the ratio 
of splittings will shift. From potential models one can argue that the IS-IP splitting will suffer 
a larger change than the 1S-2S splitting, and that the ratio, R, defined above will decrease as 
one goes from the quenched to unquenched world. The other two major sources of systematic 
error in our calculation are fitting errors and the finite lattice spacing and relativistic errors in 
the NRQCD propagators. Our experience with fitting indicates that the results tend to move 
around by about a standard deviation when different methods are used, in other words the 
systematic fitting errors are probably about the same size as the statistical errors, while we 
expect systematic errors in the propagators to contribute at the 5 MeV level to the spectrum. 

In order to present dimensionful results we need to define one global ( or "average") a~^. 
We do so by calculating a bootstrap average based on the 1S-2S and IS-IP splittings. The 
results are shown in Table 0. There we give a~^'s, separately for the two splittings plus the 
bootstrap values for averages and differences. We show results for two different ensembles for 
the P-states, the first obtained using smeared-local fits and the second based on matrix fits ( 
see Table | for spread in results due to fitting procedure dependence)). We use the bootstrap 
average plus its associated error to define a global a~^: 

a-^ = 2A{l)GeV 

Using the central value = 2.4GeV, one can convert Table ^ into physical numbers. In 
Table |Tl] we do that for aMjj = 1.71. and compare with experiment. The agreement between 



simulation results and experiment is excellent, with most entries agreeing within la. It is 
hoped that once unquenched calculations are completed, lattice QCD results will follow ex- 
periment even more closely. Our experience with the present quenched calculations tells us 
that unquenching effects are small (this observation also follows from potential model calcu- 
lations) and that high statistics data are required to see them. The entries in Table |Tl| that 
do not have a corresponding experimental number attached are predictions of the theory. We 
predict T D-states with center of mass at (10.20 ± 0.07 ± 0.03)GeV , where the second error 
comes from uncertainties in a~^. Our current simulations give an rjb state at (9.431 ± 0.005 
± 0.001)GeV, where now the first (and dominant) error is due to higher order relativistic 
and finite a corrections not included in our action and the second error is due to a~^. The 
rjh energy is likely to change (decrease) when one goes beyond the quenched approximation. 
Given the small statistical errors on the rjf, - T(IS') splitting this is one place where effects of 
quenching could be observable. 

The spectrum results are shown in Figures ^and |^ of section 1. The P-state fine structure 
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from Table ^ 

a'^{lS-2S) 


a'\lS - IP) 


Bootstrap 
a-\lS-2S) 


Results 
a-^{lS-lP) 


Average a ^ 




1.71 


2.33(11) 


2.55(12) 


2.34(8) 


2.49(4) 


2.42(6) 


0.15(7) 










2.58(14) 


2.46(9) 


0.23(13) 


1.8 


2.35(11) 


2.49(17) 


2.32(10) 


2.49(11) 


2.41(9) 


0.17(10) 










2.50(15) 


2.41(10) 


0.18(15) 


2.0 


2.39(11) 


2.51(15) 


2.39(12) 


2.50(10) 


2.45(10) 


0.12(10) 










2.52(14) 


2.45(11) 


0.13(14) 



Table 10: a'^ from 1S-2S and IS-IP splittings 





Simulation Results [GeV] 


Experiment [GeV] 


2^Si - 




0.580(26) 


0.563 


??Sr - 




0.88(12) 


0.895 


I'Pi - 




0.414(22) 


0.440 


2iPi - 




0.86(7) 


0.800 


l'D2 - 




0.74(7) 




's, - 




0.0295(5) 






'Po 


0.048(10) 


0.053 


'P2 - 
'PcM - 


'Pi 
-'Pi 


0.019(5) 
0.005(1) 


0.021 



Table 11: NRQCD spectrum results and comparison with experiment. We use a = 2.4GeV 
and aM^ = 1.71. Systematic errors due to quenching, higher order terms in the NRQCD 
action and uncertainty in the scale are not included. 

is measured relative to the center of mass of the triplet states. We first create bootstrap 
ensembles of individual P-state energies. We then calculate bootstrap estimates for the center 
of mass and for energies relative to the center of mass. In Fig. Q the energy of the ^Pi-state 
(the hb) lies slightly below the triplet center of mass. We note, however, that within our ~5 
MeV systematic errors the two levels are consistent with each other. 

One quantity of phenomenological interest is the mesonic wave function at the origin (i.e. 
zero separation between quark and anti-quark), which for the T{mS) is 

^„(0) = <m\J2i^K^)a,xH^)\0>, 

= < m\loc >zsiz (22) 
where |0 > is the fock space vacuum, < m| is the quantum state of an T{mS) polarized in the 
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z direction, and \loc >-iS\z is the state created by our loc smearing function in the ^5*1^ channel. 
Examining Eq. (p!8D, we see that (in the hmit of an infinite number of exponentials in our fit 
ansatz) < m\loc >3s-^z measured in units of a~^/^ is just the fit parameter a{loc, m), which in 
turn is equal to ^b{loc, m). We extract b{loc,m) by fitting simultaneously to G3Sj^{loc,loc;t) 
Gs^j (1, /oc; t), 6*35^(2, Zoc; t) and (^35^ (3, /oc; t). These are fits of the same form as Eq. ([TP|), 
but including Ugc = loc. In addition to obtaining a^^^ < m\loc > directly from our fits, it can 
be obtained indirectly from our separate fits to a and b, 

^3/2 ^^^0) ^ ^ ^3/, <n^,\m><m\loc>^ 

a{nsc,Tn) < nsclm > 

for any initial smearing function n^c. In Table ^we have tabulated measurements of a^^^ \I'm(0), 
from direct fits to b{loc, m) and from ratios of b{nsc, to) and a{nsc-, ^) taken from Table | and 
Table ||. We also present results for 'qh- For ground states, values obtained using different 
initial smearing functions are consistent. In general, we believe the diagonal entries in Ta- 
ble |T2| with Use = fn are the most reliable. Good smearing functions imply that off-diagonal 
amplitudes are severely suppressed. Signal to noise for the off-diagonal amplitudes a{n, m) 
and b{n, m) and for their ratios is not as good as for the diagonal amplitudes. The "direct fit" 
estimates for ^1^(0) suffer from the fact that 6(/oc, m) comes from fitting the local-local meson 
correlation function. This correlation has many more exponentials contributing out to large 
t-values than is the case for smeared-local correlations (this is the reason why the local-local 
correlation was not used in any of our fits to extract energies). Hence, it has proven difficult to 
extract the amplitude for any given excited state accurately using the local-local correlation. 
The indirect method of Eq. (p3D gets around this problem and uses only smeared-local and 
smeared-smeared correlations. 

Care should be used in any attempt to relate these values of the wave function to predictions 
of physical processes, such as the leptonic width of the T. A complete calculation of such 
processes will involve large (10 - 20% or more) corrections in relating the lattice current 
appearing in Eq. (|22D to continuum currents. These corrections arise both from relativistic 
(coming from the small components of the heavy Dirac spinors) and renormalization effects. In 
addition, |\I'(0)p scales like a~^; meaning that lattice spacing errors in the dimensionful value 
of |\I'(0)p will be of order 13%. Quenching effects are likely to induce an additional systematic 
error on the simulation results, which we have not attempted to estimate quantitatively. If 
one nonetheless goes ahead and combines the numbers in Table ^ with = 2.4(l)GeV and 
with the standard leading order Van Royen-Weisskopf formula for a vector meson leptonic 
width []T^, one finds for the T(IS'), Fee = (1.09 ± 0.03 ± 0.14)keV. The experimental value is 
T^fP'\T) = (1.34 ± 0.04)keV. In the simulation result we only quote two errors corresponding 
respectively to statistical errors and a^^ uncertainties. Quenching errors and relativistic and 
matching corrections (corrections that are expected to dominate) have not been included. 
Unquenching should increase the theoretical estimate for Fee by enhancing |\E'(0)p. Using 



the diagonal entries in Table |T2|, one finds for the T 2S and 3S, Fee(25') = (0.66 ± 0.10 ± 
0.09)keV and Fee(3S') = (0.75 ± 0.37 ± 0.10)keV. The experimental numbers are respectively 
(using the branching ratios), T'ffP-\2S) = (0.56 ± 0.09)keV and F(fP-)(3^) = (0.44 ± 

0.04)keV. When one compares ^E'(O) given in Table |1^ for the ground and excited states one 
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\mb) 
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r\ 


Q 
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.153(2) 


.154(2) 


.162(12) 


.17(2) 


T(2S) 


.147(10) 


.16(3) 


.127(10) 




T(35) 




.32(10) 




.14(4) 




.161(2) 


.162(2) 


.184(26) 






.147(9) 


.17(4) 


.132(10) 




r]b{3S) 








.13(2) 



Table 12: Mesonic wave function at the origin, a^^^ ^(0). A dashed hne means that no signal 
could be extracted. 

finds that the ratios ^I^25(0)/^I^i5(0) and \[^35(0)/\E'i5(0) are larger than one would expect from 
lattice potential calculations. Our excited state leptonic widths lie above the experimental 
numbers ( although consistent within large errors ), whereas for T(IS') the simulation result lies 
below experiment. Once corrections to the present calculation are included ( e.g. relativistic, 
quantum loop, lattice-continuum matching, unquenching) , it will be important to monitor and 
understand how and whether these discrepancies go away. At the same time, one also needs 
to reduce the statistical errors, particularly for the excited states. We hope to be reporting 
on more accurate results in the future. 



5 T Mass and Lorentz Invariance 

As we have shown, NRQCD simulations gives accurate results for the splittings between 
bb states. These simulations can also be used to compute the full mass of the T. We use 
the T mass to determine the correct quark mass for the simulation. We have investigated two 
different methods for computing the T mass. 

The first method is to add twice the renormalized mass of the quark to the nonrelativistic 
energy i?NR(Tr) obtained from the simulation: 

Mr = 2(Z^M0 - Eo) + ^^^(T). (24) 

Here and Eq are renormalizations that are computed using perturbation theory p^] . 
Our results are shown in Table 0. The uncertainties in this procedure are due to uncalculated 
0(a^) corrections in the perturbative expansions of the renormalizations, and to uncertainties 
in a~^. 

The second method is to compute the nonrelativistic energy for T's with nonzero three 
momenta. We did this for momenta p = (0,0,1), (0,1,1), (1,1,1) and (0,0,2) in units of 47r/16a, 
and fit the resulting energies to two different parameterizations: 

Er(p) - ENRiT) = Ci^^ (25) 

tip; nrk ) ^ 1 3 V ; 
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a Enpi 


7 


a Eq 




Mr (GeV) 


1.71 
1.8 
2.0 


0.453(1) 
0.451(1) 
0.444(1) 


1.20(4) 
1.18(4) 
1.16(3) 


0.32(6) 
0.31(6) 
0.30(6) 


3.92(18) 
4.08(18) 
4.48(17) 


9.4(6) 
9.8(6) 
10.8(6) 



Table 13: The T mass as determined from the T's nomelativistic rest energy in NRQCD. The 
last column assumes = 2.4(1). 





aMkin 
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Mkin (GeV) 


1.71 


3.94(3) 
3.94(3) 


1.0(3) 


0.20(37) 


0.16 
0.19 


9.5(4) 


1.8 


4.09(3) 
4.09(3) 


1.1(3) 


0.20(35) 


0.60 
0.62 


9.8(4) 


2.0 


4.48(4) 
4.48(4) 


1.2(3) 


0.24(36) 


0.64 
0.67 


10.8(4) 


1.8 

{5H = 0) 


3.82(3) 
3.83(3) 


0.23(32) 


-0.45(41) 


0.06 
0.006 


9.2(4) 


1.8 

{uo = 1) 


4.96(3) 
4.97(3) 


1.4(3) 


0.60(35) 


0.11 
0.18 





Table 14: The T mass determined from its kinetic mass Mkm using two different fits. The 
parameters Ci and C2 are explained in the text. The last column assumes = 2.4(1). 



and 

In a Lorentz invariant theory, the kinetic mass Mkin is equal to the rest mass of the upsilon; 
in our simulations this should be true up to corrections of order v^, except in the 6H = 
theory where the errors should be order f^. The (p^)^ terms in our fits test for higher- 
order relativistic effects; we expect Ci = 1 up to corrections of order (9(f^,a^), except in the 
6H = theory where this parameter should almost vanish. The pi term in the second fit tests 
for contributions that are not rotationally invariant; we expect C2 = up to corrections of 
0{v'^, a^). These expectations are confirmed by our results which are shown in Table |1^. 

Our two methods give identical T masses in lattice units to better than 1%. Assuming 
a^^ = 2.4(1), our results indicate that the correct bare mass for a b quark is M° = 1.7(l)/a; 
among our simulations, the set with aM° = 1.71 is the best. 
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The close agreement between our two determinations is striking confirmation of the partial 
restoration of Lorentz invariance due to the correction terms 5H that we include in the quark 
evolution equation. The kinetic mass is quite sensitive to these correction terms; it is too 
small by almost 10% when we set 5H — 0. This is what we expect since Mkin equals twice 
the quark mass in a nonrelativistic theory, while in a relativistic theory the binding energy 
contributes as well. Note that tadpole improvement is also essential; without it, Mkin is too 
large by almost 25%. 

Each of the two methods we use to determine the T mass has its strengths and weaknesses. 
Using Mkin is attractive because there is no need for perturbative calculations, but it is only 
accurate if the 5H correction terms are included in the quark propagators. The rest energy of 
the meson is quite insensitive to 5H , but can only be used when combined with perturbative 
results for the quark mass renormalizations. Ideally one uses both methods and compares 
results, as we have done here. 



6 Summary 

We have investigated the T system using nonrelativistic lattice QCD (NRQCD) and find that 
one can successfully reproduce the general features of the known T spectrum, including P- 
state fine structure and several excited states. We also make predictions for D-states and for 
776. By reducing the statistical errors in the splittings by a modest amount, one should, in the 
future, be able to study effects of quenching in a quantitative way. 

We worked with an action correct through 0{Mi,v^) and whose link variables had been 
divided by mq, the fourth root of the plaquette value, i.e. with a tadpole-improved action. 
Because of tadpole-improvement we were able to get good results (e.g. the correct P-state 
fine structure splittings) with tree-level coefficients for the correction terms in the NRQCD 
action. Our program should be viewed as a demonstration of a successful implementation of 
perturbative improvement of a lattice action. Our results underscore the utility of tree-level 
improvement and the crucial importance of tadpole-improvement. 

We have put considerable effort into developing fitting procedures that allow us to get 
not just ground state energies and energy splittings between levels with different quantum 
numbers, but also (radially) excited energy levels. We find that it is important to fit several 
correlations simultaneously, and do multi-exponential multi-correlation fits. 

Our calculations have systematic errors at the 5MeV level coming from higher order rel- 
ativistic and finite lattice spacing effects. The largest source of uncertainty comes, however, 
from the quenched approximation. Calculations involving dynamical configurations are un- 
derway, and we will to be reporting on them soon. Quenched cc spectrum results will also 
appear shortly. 
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